Vibrations and diverging length scales near the unjamming transition 
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We numerically study the vibrations of jammed packings of particles interacting with finite-range, 
repulsive potentials at zero temperature. As the packing fraction <f> is lowered towards the onset of 
unjamming at (j> c , the density of vibrational states approaches a non-zero value in the limit of zero 
frequency. For <j> > <f> c , there is a crossover frequency, uj* below which the density of states drops 
towards zero. This crossover frequency obeys power-law scaling with <f> — <f> c . Characteristic length 
scales, determined from the dominant wavevector contributing to the eigenmode at uj* , diverge as 
power-laws at the unjamming transition. 
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The jamming/ unjamming transition for zero- 
temperature sphere packings has a mixed first- 
order/second-order character Q, 0]. This transition 
occurs for particles with finite-range and purely repul- 
sive interactions. As a system unjams with decreasing 
packing fraction, the number of interacting neighbors 
per particle drops discontinuously to zero. Despite 
this characteristic signature of first-order behavior, 
power-law scaling is also observed for other quantities 
This raises the question of whether there is a 
diverging length scale associated with the loss of rigidity. 
Simulations suggest that a diverging length scale exists 
on the low-density side of the transition 0, but 
there has been no demonstration of similar behavior 
in the jammed phase. Because the jamming transition 
may control the glass transition at higher densities 
and temperatures Q, an observation of a diverging 
correlation length should shed light on the nature of 
the glass transition and properties of glasses in general. 
Here we show that a diverging length can be extracted 
from the vibrational properties of the jammed phase. 

In most crystalline or amorphous solids, vibrations at 
low frequency, w, are expected to be long- wavelength, 
acoustic plane waves. From this assumption one obtains 
the asymptotic low-frequency Debye density of vibra- 
tional states: T>(uj) oc uj d ~ 1 where D is the dimension 
of space. A feature of glassy systems is that there is 
an excess in the density of low-frequency vibrations Q. 
There has been an extensive body of simulation work 
on phonon spectra in quenched glasses dating back to 
the 1970's _5|. In particular, studies of the Lennard- 
Jones glass found an increasing number of low-frequency 
modes as the system was diluted in an ad-hoc fashion [6j. 
Perhaps the most dramatic demonstration of this excess 



density of low-frequency modes appears in a system at 
zero temperature with a well-defined jamming transition 
at packing fraction 4> c 0]. For systems with attractions, 
the jamming transition is usually inaccessible because of 
the vapor- liquid phase transition Q. Therefore, a sys- 
tem with soft, purely-repulsive, finite-range interactions 
has the advantage of allowing a systematic study of the 
properties near its well-defined jamming transition. 

In the simulations reported here we have studied 
monodisperse spheres of diameter er interacting through 
a finite-range, purely repulsive, harmonic potential: 



V(r) = { 







r < a 
r > a. 



Here r is the center-to-center separation between two 
particles. The units of length and time are a and 
(md 2 /Vo) 1 / 2 , respectively, for particles of mass m. Our 
three dimensional (3D) systems consist of N = 1024, 
4096, and 10000 spheres in cubic simulation cells with 
periodic boundary conditions. We employ conjugate- 
gradient energy minimization in order to obtain T = 
configurations. We have also studied bidisperse, har- 
monic discs in 2D, as well as Hertzian contact poten- 
tials in 3D. Our conclusions are the same for all three 
systems. 

In Fig. ^a) we show a log-log plot of the density of 
states, T>(uj) as a function of uj covering eight decades in 
(4> — <f) c ) for one system size. For the smallest value of 
(<j)—<f) c ) studied, the low-frequency behavior approaches a 
flat spectrum with an uj — intercept of Vq = T>(uj — > 0). 
Thus, close to the jamming transition there is no longer 
any vestige of Debye behavior. As the system is com- 
pressed above threshold, the curves depart from this 
plateau behavior at a frequency uj* , which increases with 
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(4> — 4> c ). Below to*, T>(uj) — * as u — > 0. Experiments 
on vitreous silica || and simulations of models of glasses 
with three-body interactions |9( have observed a similar 
trend with decreasing pressure. Figure IHb) is the full 
spectrum divided by lu 2 for all system sizes, at several 
values of (<f> — </> c ). The peak position in T>(lu)/lu 2 shifts 
to lower frequencies and the peak height increases as 4> c 
is approached, analogous to the way some glasses behave 
as the glass transition temperature is approached |ld| . 
Thus, our zero-temperature system captures features of- 
ten associated with the Boson peak in real glasses [llj . 
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FIG. 1: (a) Density of vibrational states T>(u>), for N = 1024 
monodisperse spheres in 3D, interacting via harmonic re- 
pulsions. The curves from right to left are for, <f> — <f> c — 
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4>-<i>c = lx 10~ 6 , lxlO -4 , lxl0 _s , 1x10 
using data from N = 1024, 4096, and 10000. 



(b) V(uj)/lo 2 at 
~ 2 , and 1 x 10~\ 



Figure[21shows the crossover frequency, lu* , versus (</> — 
4> c ). We determine lu* by finding where "D(ui) for a given 
{4> — <j) c ) departs from T>(uS) for the next smallest value 
of ((f) — (f> c ). Over more than 4 decades in (<f> — cf> c ), lu* 
vanishes as a power law: 



(1) 



with Q = 0.48 ±0.03. The scaling is robust, independent 
of the precise manner in which we determine to* . It is the 
same (with different prefactors) when u* is defined as the 



value of a; at which T>(u>) reaches Dq, 0.9 Do or 0.5£>0j as 
well as the frequency where T>(u>) /lu 2 in Fig. Cfb) levels 
off with decreasing w. 
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FIG. 2: The crossover frequency ui* . The line is a fit to Eq.0 
with fi = 0.48. 

Given lu*, we can define a corresponding length scale, 
£r, that diverges as <p approaches (f>J. In order to ex- 
tract we examine the spatial variation of the eigen- 
mode corresponding to the frequency lu*. We take the 
Fourier transform of the appropriate component of the 
polarization vector, P,-(o;*) of each particle i, which, for 
transverse waves is, 



Mk,Lu*) 



5^fcAPi(w*)exp(*kTi) 



where ( } denotes an average over configurations and over 
all wave- vectors with the same magnitude k. (The lon- 
gitudinal component, /z,(fc, to), not shown ^3], is the dy- 
namical structure factor accessible from inelastic neutron 
scattering.) Figure^a) shows these functions at the val- 
ues of lu* determined for different compressions, (<j) — <j) c ). 
The transverse components exhibit well-defined peaks at 
k* at small wavevectors (there is a system size cut-off at 
k m in — 27r/L, where L is the size of the simulation box). 
Thus, £t = 27r/fc* is the dominant transverse length scale 
for that mode. Figure EHb) shows £t as a function of 
(4> — 4> c )- The solid line corresponds to a power-law fit: 



£t oc 



(2) 



with v T = 0.24 ± 0.03. The exponents v T and fi (Eq.QJ 
can be related to each other via a simple scaling argu- 
ment using k* = lu* /op (</>), where ct(</>) is the velocity of 
transverse sound. The shear modulus Goo vanishes with 
an exponent of 7 = 0.48 ± 0.05 |2|, implying that ct(4>) 
vanishes with an exponent of 0.24. Given = 0.48, we 
would therefore expect vt = H — 7/2 w 0.24, in agree- 
ment with the results of Fig. Otb). 

A similar argument can be constructed for a longitudi- 
nal length £l based on the peak of /r,(fe,o;*). Although 
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FIG. 3: (a) Transverse spectral functions fr{k, u)*) for <f) — 
4> c = lxlO" 5 , 5xl(T 5 , lxlO" 4 , 5xl0" 4 , lxlO" 3 , 5xl(T 3 , 
with decreasing amplitude respectively, (b) The correlation 
length £t = 2n/k* obtained from the wave- vector k* for the 
position of the peak in fx at frequency u*. 

it is difficult to extract the peak value from f^(k,u>*) be- 
cause the peaks occur at very low values of k, the anal- 
ogous scaling relation based on the longitudinal sound 
speed (and hence the bulk modulus, which is indepen- 
dent of 4> — (j) c 0) predicts 

v L « 0.48. (3) 

We have also simulated a system with Hertzian po- 
tentials: V {1 - r/cr) 5 / 2 , for r < a, and zero otherwise. 
In this case, we find that to* vanishes with an exponent 
SI « 0.78 ±0.03, which is different from that obtained for 
the harmonic case. However, by calculating the peak of 
fr(k,uj*) we obtain a length scale £t that diverges with 
vt ~ 0.23, just as in the harmonic case. This is again 
consistent with the scaling argument based on the speed 
of sound, since the shear modulus vanishes with an ex- 
ponent 7 w 0.95 for the Hertzian case |2(. Thus, even 
though the scaling of u>* is different for the Hertzian and 
harmonic potentials, we have the same scaling for £t. 

The strong departure from the low-frequency Debye 
density of states suggests that the eigenmodes are poorly 
characterized by simple plane waves ^3|- We illustrate 



this point in Fig. ^ where we show the lowest frequency 
modes for 2D harmonic systems at the two extreme val- 
ues of O - (f> c ) = 1 x I0" 1 and 1 x I0~ 8 . We show 2D 
results here for ease of visualization. These correspond to 
modes below and above u>* , respectively, for those pack- 
ing fractions. In the compressed system at u> < oj* the 
eigenmode retains identifiable correlated plane-wave-like 
character, consistent with the more Debye-like behavior 
of T>{u>) as lo — > 0. In contrast, when lo > cu* as in the 
system closest to <f> c , all plane- wave character is lost. This 
picture suggests that £t represents the length scale above 
which one can average over disorder: wavevectors of mag- 
nitude k > 27t/£t do not have long enough wavelengths 
for effective averaging to take place. This is related to the 
criterion that characterizes the elastic-granular crossover 
in a Lennard- Jones glass [lij . 




FIG. 4: Lowest frequency eigenmodes for N = 1024 bidisperse 
discs in 2D, at two extreme values of <j> — <j) c : 1 x 1CP 1 (left), 
and 1 x 10 -8 (right). The dots represent the centers of the 
particle and the black lines their polarization vectors. 

Our results for uj* and £r motivated two recent the- 
oretical papers that treat different aspects of the zero- 
temperature jamming transition. The first paper, by 
Wyart et al. ^^|, shows that the constant density of 
states is a direct consequence of the system being at what 
is called an isostatic point, where the number of con- 
strai nts p recisely equals the number of degrees of freedom 
[1HG3- This theor y Predicts CI = 1/2 and Q = 3/4 for 
the harmonic and Hertzian potentials, respectively, in ex- 
cellent agreement with our simulations. It also predicts 
a new length scale, £ oc {(f> — § c )~ x l 2 . This is not the 
same scaling as we have found for but it is consistent 
with the scaling of £l predicted by the speed of sound 
argument (Eq. |3J). The length scale defined in Ref. ^| 
results from a competition between bulk mechanical sta- 
bility and boundary-induced mode softening, and has not 
yet been explored by simulation. 

A second theory by Schwarz et al. 17] makes an anal- 
ogy between jamming at cj> c and the onset of fc-core per- 
colation. The fc-fold coordination required for each site 
corresponds to the D + 1-fold coordination required for 
a particle to be locally stable in D dimensions. For the 
Bethe lattice, Schwarz et al. [l7j find exponents for the 
number of overlapping neighbors per particle and the sin- 
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gular part of the shear modulus that are in good agree- 
ment with our simulations [lj. Moreover, they find a cor- 
relation length exponent of v — 1/4, in agreement with 
Eq. |21 Interestingly, the Bethe lattice calculations 0] 
also suggest the presence of another length scale with an 
exponent of v# = 1/2, in possible agreement with pre- 
dictions of Wyart et al. |l5| and our argument for £r, . 

O'Hern et al. Q, Q found another correlation- length 
exponent by measuring the shift in the peak of the dis- 
tribution of jamming thresholds as the system size in- 
creases 0,0. This exponent was measured to be v « 0.7 
in both 2 and 3 dimensions. Finally, yet another di- 
verging length scale has been measured in simulations by 
Drocco et al. Q as the jamming threshold is approached 
from below. They measured the distance over which a 
slowly-moving particle disturbs the surrounding packing 
and found V- = 0.6 — 0.7, in good agreement with the 
finite-size scaling results of O'Hern et al. 0,0- 

There are at least three length scales that are impor- 
tant in describing the jamming/ unjamming transition: 

(i) The intcrparticle overlap distance, which goes to zero 
at (fi c with a power-law exponent of approximately 1.0 Q. 

(ii) the length scale determined from finite-size scaling 
y], as well as the length Q calculated as the transition 
is approached from low densities. These diverge with an 
exponent of roughly 0.7. (iii) the diverging length scale 
presented here on the high-density side of the transition, 
determined from transverse vibrations; this diverges with 
an exponent of 0.24. Theories and a scaling argu- 
ment based on our data (Eq. |3J) suggest a fourth length 
scale with exponent 0.5. It is still not clear how all these 
different length scales are tied together. That different 
exponents are observed on different sides of the transition 
suggests that the jammed phase may always be separated 
from the unjammed one by a singularity. 

Disordered or glassy systems characteristically show 
an excess in the number of low-energy excitations. In 
particular, there is a Boson peak that can increase in 
magnitude and shift to lower frequency as the glass tran- 
sition is approached 0, ^| . However, nowhere are these 
two effects more clearly demonstrated than at the on- 
set of unjamming. There T>(u>) approaches a constant as 
lu — > instead of dropping to zero, implying, as shown in 
Fig.^b), that the Boson peak diverges and the peak po- 
sition, u)* , shifts all the way to zero frequency as the rigid- 
ity disappears. This is in accord with random-network 
theories |l8j and numerical studies imposing disordered 
force constants on a lattice • The behavior we observe 
is clearly associated with the critical nature of the transi- 
tion and thus provides a natural way of describing these 
effects in terms of a diverging length scale. Our studies 
also indicate that the enhancement of the low-frequency 
modes is purely a geometrical phenomenon. 

In conclusion, we have studied the properties of the 
jamming/unjamming transition at zero temperature. 
Despite the discontinuity in the number of interacting 



neighbors, we find that the loss of rigidity is characterized 
by a diverging length scale. Moreover, we note that the 
jump in the number of interacting neighbors is universal, 
and given by the isostatic condition. This indicates that 
unjamming may be described as a second-order transi- 
tion with a universal jump in the order parameter. 
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